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Abstract 

The  response  of  vertical  arrays  at  single  frequencies  (CW)  and  for  homogeneous  media 
is  well  known.  This  paper  addresses  the  issues  of  frequency  dependence  and  sound 
velocity  gradients  for  the  vertical  array  response  in  a  deep  ocean.  I  have  modified 
the  synthetic  seismogram  code  of  Neil  Frazer,  Subhashis  Mallick  and  Dennis  Lindwall 
to  address  this  problem.  The  code  uses  a  rearrangement  of  the  Kennett  reflectivity 
algorithm  (Kennett,  1974,  1983)  which  computes  the  geoacoustic  response  for  depth 
dependent  media  and  pulse  sources  by  the  wave  number  integration  method.  The 
pneraHzed  Filon  method  is  appHed  to  the  slowness  integral  for  an  additional  increase 
in  speed  (Frazer  and  Gettrust,  1984;  Filon,  1928).  The  original  code  computes  the  re¬ 
sponse  of  a  single  source  at  a  specified  depth.  The  new  code  has  several  improvements 
over  the  previous  one.  First,  it  is  a  much  simplified  code  addressing  only  acoustic 
interaction.  The  total  length  is  about  half  the  length  of  the  original  code.  Secondly, 
the  code  can  compute  the  response  of  a  vertical  array  of  point  sources.  By  changing 
the  phase  delay  between  the  sources,  we  can  steer  the  beam  to  the  places  of  most 
interest.  Thirdly,  the  code  reduces  considerably  numerical  noise  at  large  offsets.  The 
original  work  has  numerical  noise  beyond  about  30  km  offset  at  50  Hz  which  bmits 
the  application  of  reflectivity  modeling  in  long  range  problems.  The  improvement 
comes  with  the  optimization  of  the  program,  both  in  the  speed  and  program  struc¬ 
ture.  The  improved  algorithm  can  be  used  to  get  the  far  offset  response  (up  to  150 
km)  of  a  vertical  array  in  the  deep  ocean  at  frequencies  up  to  at  least  250  Hz.  The 
modeling  results  are  compared  to  analytical  and  benchmark  solutions.  The  modified 
reflectivity  code  can  be  applied  to  the  study  of  pulsed- vertical  array  sources  such  as 
were  deployed  on  the  ARSRP  (Acoustic  Reverberation  Special  Research  Program) 
acoustic  cruises. 
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Chapter  1 


Introduction 


The  ARSRP  (Acoustic  Reverberation  Special  Research  Program)  Reconnaissance  Ex¬ 
periment  was  conducted  from  25  July  to  19  August  1991,  and  in  1993  the  ARSRP 
acquired  detailed  geological  and  acoustic  backscatter  data  from  three  sites  in  the 
ARSRP  corridor  in  the  western  North  Atlantic.  A  vertical  line  array  (VLA)  of  10 
coherent  sources  was  used,  and  its  beam  angle  can  be  steered  (EUiot,  1991).  Spe¬ 
cific  long  term  objectives  of  these  efforts  are;  1)  “to  characterize  the  variations  in 
bottom  topography  and  sub-bottom  properties  that  control  the  scattering  of  low  fre¬ 
quency  acoustic  waves”,  2)  “to  develop  theoretical  and  numerical  techniques  capable 
of  predicting  the  low  frequency  acoustic  wavefield  scattered  from  geologically  realistic 
models  of  the  bottom/sub-bottom  environment”,  and  3)  “to  isolate  from  these  scat¬ 
tering  models  the  physical  mechanisms  which  dominate  the  long-range  reverberation 
from  the  seafloor.” 

The  challenge  for  the  ARSRP  is  to  analyze  the  acoustic  and  geological  data  and 
to  demonstrate  a  predictive  modeling  capability  for  low  angle  seafloor  backscatter. 
In  order  to  achieve  such  a  goal,  it  is  important  to  model  the  acoustic  responses  on 
the  seafloor  for  the  VLA.  There  are  several  interesting  objectives  in  the  numerical 
modeling:  1)  what  is  the  effect  of  the  free  surface  on  the  sea  floor  response,  2)  what  is 
the  effect  of  the  ocean  velocity  gradient  on  the  response,  3)  how  do  the  above  factors 
affect  CW  and  pulse  sources,  etc. 

Ray  theory  and  parabolic  equation  methods  (Smith  and  Tappert,  1993)  are  com- 
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mon  in  ocean  acoustics  to  compute  the  acoustic  field.  They  are  compared  in  Tolstoy 
et  al.  (1985).  Ozliier  (1992)  studied  the  refraction  effects  on  vertical  line  array 
beamforming  applying  a  simple  ray  theory  method.  She  studied  the  responses  from 
10  omnidirectional  point  sources  with  linear  phase  tapering  equivalent  to  a  steering 
of  6  deg  depression.  There  are  big  inaccuracies  involved  in  the  results  after  30  km 
horizontal  offset. 

In  order  to  study  the  vertical  array  interference  problem  more  completely  and 
to  get  a  more  accurate  picture  of  the  interference  response  in  the  deep  ocean  for  a 
wide  variety  of  outgoing  beam  angles  from  the  vertical  array,  we  use  the  reflectivity 
modeling  method.  The  reflectivity  method  has  been  widely  used  to  compute  synthetic 
seismograms  in  layered  media.  In  fact,  it  has  contributed  to  a  better  understanding 
of  the  earth  s  structure,  both  on  the  continents  and  beneath  the  ocean  (for  example, 
Braille  and  Smith  (1975),  Spudich  and  Orcutt  (1980a,  b),  and  Kempner  and  Gettrust 
(1982a,b)). 

Its  main  advantage  is  its  capacity  to  compute  a  total  solution  of  the  wave  field 
for  a  given  model.  A  matrix  method  is  generally  used  to  compute  the  response 
of  the  model  in  frequency- wavenumber  space.  This  includes  contributions  from  aU 
possible  generalized  rays  within  the  reflecting  zone  (Kennett  1974,  1983;  Kind,  1976). 
The  original  reflectivity  version  of  Fuchs  and  MiiUer  (1971)  required  the  source  and 
receiver  to  be  above  the  reflecting  zone,  but  the  method  was  subsequently  modified 
by  Stephen  (1977)  to  accommodate  a  receiver  buried  within  the  reflecting  zone.  In 
practice,  thei^  are  some  disadvantages  of  the  reflectivity  method.  The  main  one  is 
the  usually  long  computation  time  required  in  the  modeling.  So,  even  though  there 
has  been  extensive  study  of  the  theory  of  the  reflectivity  method,  adequate  care  must 
be  taken  to  implement  the  theory  and  also  to  balance  the  speed  and  accuracy.  This 
is  especially  true  for  our  vertical  array  problem  which  has  multiple  sources  and  which 
is  required  to  compute  long  range  responses  in  the  deep  ocean. 
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Chapter  2 


Background 

2.1  Reflectivity  Method 

The  reflectivity  method  has  long  been  used  by  seismologists  for  modeling  both  land 
and  marine  reflection  and  refraction  data.  The  method,  originally  proposed  by  Fuchs 
and  Muller(1971),  was  extensively  modified  by  Kennett  (1974).  The  calculations  are 
done  essentially  in  two  steps: 

•  A  reflectivity  function  iZ(a;,p)  is  calculated  in  the  frequency-ray  parameter 
(wavenumber)  domain.  This  is  performed  by  layer  iteration  starting  from  the 
free  surface  down  to  the  deepest  interface.  In  each  step  of  the  iteration,  all 
orders  of  multiple-bounce  paths  in  the  layer  are  included.  At  the  final  step, 
one  obtains  a  reflectivity  function  that  includes  aU  possible  ray  paths  from  the 
source  to  the  receiver. 

•  The  second  step  involves  numerical  evaluation  of  a  double  integral  of  the  form 

da;  exp(— tc<;i)  X 

where  cr  ^  iu)x  and  g{jp)  ==  p. 

The  integration  over  frequency  a;  is  usually  carried  out  by  a  fast  Fourier  transform 
(FFT).  A  complex  frequency  with  a  constant  imaginary  part  is  used  in  the  integral 


j  <^P/(w,P)exp[<T£7(p)]  (2.1) 


uix,t)  =  - 
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to  attenuate  the  wraparound  caused  by  the  use  of  the  FFT. 

The  integrand  of  the  integral  over  the  ray  parameter  p  is  highly  oscillatory  and  use 
of  the  trapezoidal  rule  would  require  a  very  small  step  size  in  p.  Use  of  a  generalized 
Filon  method  (GFM)  (  Frazer,  1978;  Frazer  and  Gettrust,  1984)  allows  one  to  use  a 
much  larger  step  size  in  p.  The  sampling  interval  5p  depends  both  on  frequency  and 
the  maximum  range  needed  for  the  calculation,  i.e.,  long  range  and  high-frequency 
calculations  require  a  very  small  step  size  in  p  to  avoid  spatial  aliasing. 

The  background  and  limitations  of  the  method  are  described  in  the  two  papers 
of  Mallick  and  Frazer  (1987;  1988).  We  modified  their  program  for  our  vertical  array 
problem. 

The  compressional  potential  of  the  wave  from  an  explosive  point  source  is 

MT,z,t)=^F{t--)  (2.2) 

It  Oci 

where  =  r*  +  z*.  Its  Fourier  trsuasform  can  be  written  in  integreJ  form 

_  _ 

^Jr,z,<v)  =  F(u})  I  — Jo(kr)exp(— ji/iz)dA:  (2.3) 

Jo  JVi 

where  F{uf)  is  the  Fourier  transform  of  the  excitation  function  F{t),  Jo(fer)  the  Bessel 
function  of  the  first  kind  and  order  zero,  j  the  imaginary  unit,  k  the  horizontal  wave 
number,  and 

=  (K,  -  (2-*) 

is  the  vertical  wave  number  (kai  = 

Since  we  are  mainly  interested  in  the  application  of  synthetic  seismograms  to 
explosion  seismological  studies,  we  can  replace  the  Bessel  functions  by  their  asymp¬ 
totic  approximations  for  large  atguments(Fadis  and  Muller,  1971)  ,  which  is  good 
for  kr  >  14  (Corresponding  to  a  source  with  a  frequency  of  200  Hz,  this  approxi¬ 
mation  is  good  for  ranges  down  to  16.7m.  The  synthetic  seismograms  based  on  this 
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approximation  will  be  incorrect  at  normal  incidence)  : 


[Y 

1 

f  r  *  /  \  T 

TT 

1 

1/0(2^) 

~  \  —  cos(x  - 
V  TTX 

's/2'kx 

|_exp[;(a;  -  -)]  -b  exp[-j 

'■'^-4)1 

[  (2.5) 

Y,{x) 

[Y  . 

%  \  —  sinfa:  — 
V  Tra: 

-  ] 

TT 

exp[j(x  -  -)]  -  exp[-j 

>  (2.6) 

Where  Jq  and  Yq  are  Bessel  functions  of  first  and  second  kind  respectively  and  both 
are  of  order  zero. 

The  second  exponential  term  in  the  above  equations  corresponds  to  waves  prop¬ 
agating  in  the  positive  r-direction  (away  from  the  source),  whereas  the  first  term 
describes  waves  traveling  in  the  negative  r-direction  (towards  the  source). 


2.2  The  Reflectivity  Function  in  a  Layered  Water 
Column 

Fig.  2-1  shows  the  geometry  and  notation  for  the  derivation  of  a  layered  half  space. 
The  pressure  reflection  and  transmission  coefficients  at  the  boundary  of  layers  1  and 
2  for  a  wave  incident  from  above  are: 


5?12(^i) 


P2C2  COS  0i  —  Pi  Cl  cos  02 
P2C2  cos  01  -b  Pi  Cl  cos  02 


^12(^1) 


2p2C2  cos  01 

P2C2  COS  01  +  Pi  Cl  COS  02 


with  similar  expressions  for  3f?23  and  S23.  From  Figure  2-1,  the  total  up-traveling 
signal  is  the  sum  of  an  infinite  number  of  partial  transmissions  and  reflections.  Each 
path  within  the  layer  has  a  phase  delay  2/52^2  cos  ^25  where  ^2  cos  ^2  is  the  vertical 
component  of  the  wave  number  in  the  layer.  By  letting  the  incident  signal  have  unit 
amplitude,  the  total  reflection  3fJi3  is 
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Figure  2-1;  Reflection  and  transmission  at  a  thin  layers  From  Figure  2.10.2.  in  Clay 
and  Medwin.  1976) 


—  3^12  "F  ^12^21^23  ^■^p( — 2i$2)  +  Si2^2i3?23  3i21  4i$2)  "T  ■  •  •  (2*9) 


#2  =  kih2  00^62  (2.10) 

After  3?i2,  terms  in  (  2.9)  have  the  form  of  a  geometric  series 

CO 

5  =  ^  r"'  =  (1  —  r)~^  for  r  <  1 

n=0 


3?13  =  3?, 2  +  ;?12S2l3?23eXp(-2i$2)E  [^233^21  e.xp(-2i<i2)r  (2.11) 

0 

Note  that  the  reflection  and  transmission  coefficients  at  a  single  interface  are 
frequency  independent  (i.e.  a  function  of  angle  or  ray  parameter  only).  When  the 
propagation  through  a  layer  is  considered  (i.e.  using  $2)1  the  coefficients  become 
frequency  dependent. 

We  can  reduce  equation  2.11  by  using  the  following  relations,  which  come  from 
equations  2.7  and  2.8. 
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3?12  ~  — 3?21 


(2.12) 


We  have 


^12^21  —  1  —  3?12^. 


3?12  +  3^23  exp(  — 2i$2) 


1  +  5?123?23  exp(-22$2)' 

The  transmission  through  the  layer  for  a  unit  incident  signal  is 


(2.13) 


(2.14) 


^13  “  ^12^23  6Xp(—i$2)  +  ^12^233?233?21  exp(  — 3z$2 


(2.15) 


This  is'a  geometric  series,  and  the  sum  is 


Q:  _ 


%2%3  exp(-z$2) 

1  +  §ffi23?23  exp(— 22$2y 


Both  the  3?i3  and  S513  are  oscillatory  functions  and  depend  on  $2  =  ^12^2  cos  ^2- 
They  are  functions  of  frequency  and  angle  of  incidence  for  a  given  layer. 

We  then  can  derive  the  total  reflection  and  transmission  of  n  layers,  by  repeated 
applications  of  the  single  layer  coefficient. 

As  in  Fig.  2-2  the  reflection  from  the  lower  half  space  is  Applying  equa¬ 

tion  2.14  the  reflection  coefficient  at  the  top  of  the  n  —  1  layer  is 

SfJ  _  ^(n-2)(n-l)  +  jft(n-l)n  6Xp(  — 2t$„_i) 

l-b«(„_2)(n-l)3?(n-l)nexp(-2i$„_i)'  ^ 

We  can  repeat  the  above  process  to  get  §?(n  -  3)n,  the  reflection  coefficient  for 
the  layers  beneath  the  interface,  which  is 


_  ^(n  — 3)(n— 2)  “b  — 2)n  2z$n  — 2) 

~  1  +  3?(„-3)(n-2)»(n-2)neXp(-2f$„_2)' 


(2.18) 


19 


Figure  2-2;  Reflection  from  a  layered  half  spacel  From  Figure  2.10.3.  in  Clay  and 
Medwin,  1976) 

Continuing  the  above  calculation  upward  to  the  top,  we  have 

sp  —  +  ^2nexp(— 21^2)  -Q. 

l+Si23fi2nexp(-2f$2)  ^  ’ 

By  letting  the  reflection  coefficient  3tin  represent  all  the  frequency  and  angle  de¬ 
pendence,  we  simplify  the  expression  for  the  reflection  from  a  multiple  layered  half 
space.  The  above  process  can  also  be  applied  to  get  the  transmission  coefficient  from 
a  multiple  layered  half  space.  So,  for  large  ranges,  we  can  get  the  composite  pressure 
reflection  and  transmission  coefficients  by  applying  the  above  layer  iteration  approach. 


2.3  Computation  of  the  Slowness  Integral 


In  the  reflectivity  method,  we  do  the  numerical  evaluation  of  the  slowness  integral, 
in  the  form 

yoo 

u{x,u;)  =  u^pdpJo{u;px)u{uj,p)  (2.20) 

where  Jq  denotes  the  first  kind  Bessel  function  of  order  zero.  This  integral  is  oscilla¬ 
tory  and  in  a  normal  integration  scheme  many  steps  are  necessary  in  computing  this 
when  LUX  is  large.  The  step  size  used  is  inversely  proportional  to  lux.  The  program 
uses  the  generalized  Filon  method  given  by  Frazer  (1978)  and  Frazer  and  Gettrust 
(1984),  which  requires  the  step  size  to  be  inversely  proportional  to  which 

allows  higher  u;x  response,  for  a  given  step  size  and  a  given  error  in  computation. 

To  use  the  generalized  Filon  method,  we  transform  the  integral  of  equation  2.20 
into  (see  Chapman  1978) 


^  (2.21) 

where  T  is  the  contour  of  integration  (shown  in  Fig.  2  of  Frazer  and  Gettrust  (1984)). 
The  stands  for  a  Hankel  function  of  type  1  and  order  zero.  The  details  for  the 
transformation  from  equation  2.20  to  equation  2.21  are  shown  by  Chapman  (1978). 
We  can  rewrite  equation  2.21  as 


^  /(p)e^«(p)dp, 

where 

f{p)  = 

S  —  icux 


(2.22) 


and 


gip)  =  p 
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Application  of  the  standard  trapezoidal  rule  to  the  integral  in  equation  2.22,  between 
the  limits  a  and  6,  gives  the  quadrature  formula 

I  (2.23) 

This  formula  does  not  work  well  because  it  assumes  that  is  well  approxi¬ 

mated  by  a  Hnear  function  over  the  interval  (a,  6),  while  actually  it  is  highly  oscillating. 
If,  we  assume  that  both  f{p)  and  g{p)  are  well  approximated  by  different  linear  func¬ 
tions  on  {a,b),  then  we  get  the  generalized  Filon  method  analog  of  the  trapezoidal 
rule  (Frazer,  1978): 

I'  /(p)e^«Wdp  =  -f  for  S{g)  -  0  (2.24) 

where,  Sp  denotes  p^  —  P\i  b[g)  denotes  g{p2)  —  g{pi)f  etc.  It  can  be  derived  by 
replacing  the  integrand  in  the  left-hand  side  of  2.23  by 

The  generalized  Filon  method  greatly  improves  the  quality  of  the  synthetic  result. 
It  saves  computation  time  by  as  much  as  80  percent  (MaUick  and  Frazer,  1987).  A 
straightforward  error  analysis  shows  that,  for  a  given  accuracy,  the  step  size  in  2.23 
is  proportional  to  1^1“^  whereas  the  step  size  to  2.25  is  proportional  to  |s|“^/^(Frazer 
and  Gettrust,  1984). 

The  integral  in  equation  2.25  is  simplified  by  letting  Sp  =  S{g)  and  removing  the 
term  SpfS{g): 


j  f{p)e^^^^Up  =  ■^[^(/e^®)  -  otherwise  (2.26) 

The  results  using  equation  2.26  show  great  improvement  in  terms  of  the  quality 
of  the  modeling  result.  Applying  equation  2.25  as  in  the  original  program,  there  is 
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a  big  numerical  noise  problem  beyond  about  30km  offset  which  limits  the  apphca- 
tion  of  reflectivity  modeling  in  long  range  problems.  The  optimized  program  using 
equation  2.26  reduces  significantly  the  noise  at  large  offsets.  In  the  integral  2.25,  the 
p  and  g{p)  are  very  small.  This  results  in  more  numerical  error  than  the  simpHfied 
integral  2.26. 

For  the  reflectivity  method,  we  need  to  evaluate  the  integral: 

y»00~to 

u{u},x)  =  /  kdku{u;,k)Jn{kx)  (2.27) 

Jo 

in  which  is  the  Bessel  function  of  order  n,  and  k  is  the  wavenumber.  Here  advantage 
is  taken  of  the  relation  (Olver,  1972,  formula  (9.2.19)) 

J„  =  Af„cos^„  (2.28) 

where  the  definitions  of  and  6^  are: 

M„  =  tan  “^(Fn/Jn)  (2.29) 

K  =  {Jn  +  (2.30) 

The  functions  and  On  are  available  as  polynomial  approximations  (AUen,  1954; 
Ovler,  1972,  formula  (9.4.3)  and  (9.4.6))  for  values  of  kx  greater  than  three.  For 
values  of  kx  less  than  three,  exact  values  of  Mn  and  On  could  be  computed  (Olver, 
1972,  formula  (9.2.17)). 

The  asymptotic  expansion  of  On  is  (Olver,  1972,  formula  (9.2.29)) 

On  —  kx  —  an 0{\kx\~^)  (2.31) 

where 

a„  =  \{2n  +  l)7r 
4 
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So,  the  equation  2.27  can  be  written  as: 

u{u:,x)  =  dkh{k)e'^’^  +  dkf2{k)e-'^^  (2.32) 

in  which 

h{k)  =  (2.33) 

f2{k)  =  hu{^,  A:)M„e-‘(""-*'")  (2.34) 

Equation  2.32  is  exact,  and  yet  the  functions,  fi  and  /2,  are  relatively  non-oscillatory 
because  of  2.31.  More  importantly,  each  integral  on  the  right-hand  side  of  2.32  is  of 
the  form  2.22,  and  so  they  can  be  evaluated  using  the  generalized  Filon  integration 
method. 

In  the  modeling  of  refraction  data,  x  is  usually  greater  than  four  or  five  wave¬ 
lengths,  and  then  the  following  simpler  procedure  can  be  applied.  In  equation  2.32 
replace  dj^  by  kx  —  an  and  Mn  by  (2/7rkxy^'^  (Frazer,  1988). 

The  u(a;,x)  can  now  be  evaluated  using  the  generalized  Filon  formula  2.26  with 
^(p)  =  p  and  s  =  iujx. 
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Chapter  3 


Theory  of  Vertical  Array 
Interference  for  Harmonic  Sources 


3.1  Notation  and  Example  Parameters 

In  our  study,  we  use  the  following  notation  as  shown  in  Fig.  3-1.  We  denote  W 
as  the  total  distance  between  the  top  source  and  the  bottom  source,  d  denotes  the 
neighboring  source  distance,  h  is  the  depth  of  the  top  source  from  the  free  surface 
and  (f)  is  the  grazing  or  dip  angle  of  the  ray  to  a  receiver  at  a  large  distance  from  the 
array.  Unless  otherwise  specified,  all  of  the  plotting  will  have  the  following  parameters: 
W  =  5.49A,  h  =  24.755A,  which  implies  that  the  distance  between  the  surface  to  the 
of  the  10  sources  is  27. 5A  and  the  distance  between  the  neighboring  sources 
is  0.61  A.  In  the  case  of  a  frequency  of  250jE7z  and  a  velocity  of  l.5km/s  (so  that  the 
wavelength  A  ^  6m),  d  =  3.66m  and  /i  +  ^  =  165m. 


3.2  Dipole  Interference 

3.2.1  The  Field  from  a  Single  Source 

A  sinusoidaly  excited  source  expands  and  contracts  repeatedly.  The  resulting  con¬ 
tractions  (density  increases)  and  dilatations  (density  decreases)  in  the  medium  move 
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away  from  the  source  at  the  sound  speed  c,  as  would  the  disturbance  from  an  impul¬ 
sive  source.  This  disturbance  is  called  a  continuous  wave  (CW),  and  it  comes  from  a 
CW  source.  The  distance  between  adjacent  contractions  along  the  direction  of  travel 
is  the  wavelength  A. 

Assume  that  an  omnidirectional  CW  sinusoidal  point  source  in  a  homogeneous 
medium  has  the  pressure  at  unit  distance  given  by; 


P  =  Po  *  sin(a;t) 


(3.1) 


So,  at  distance  R,  the  pressure  at  a  given  time  t  will  be: 


Pr  = 


R 


*  sina;(t  —  r) 


(3.2) 


where  r  is  the  phase  delay  due  to  the  travel  time  in  the  medium  (r  =  R/c)  and  c  is 
the  wave  velocity. 

The  field  intensity  (transmitted  power  per  unit  area),  I,  wiU  be: 

T_P^  P^ 

pc  2pcR2  (3.3) 

where  P2  the  ensemble  average  of  the  time  domain,  p  is  the  density  of 

the  medium  and  c  is  the  wave  speed  in  the  medium. 


3.2.2  The:  Field  from  Two  Sources 

Under  the  same  assumption  as  the  single  source,  two  sources  wiU  generate  an  inter¬ 
ference  pattern. 

The  pressure  at  the  receiver  R  is  given  by 


P.  =  Pi  +  P2 


where 


(3.4) 
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(3.5) 


D  -fo  .  /  ( t  \ 

Fi  =  —  *  sin  (pi  =  —  *  sin  u){t  —  Ti) 

ill  ill 

P2  =  ^  *  sin  (;^2  =  ^  *  sin  w[(i  —  T2)  +  ^o]  (3.6) 

il2  il2 

(  Ti  and  r2  are  the  wave  travel  time  from  source  1  and  source  2  to  the  receiver 
respectively,  is  the  phase  difference  of  the  sources  ) 

Then  the  field  intensity  will  be: 


/  =  =  =  +  ^  +  =  +  +  (3.7) 

pc  pc  pc  pc  pc  pc 

(ii  and  I2  are  the  intensity  due  to  single  sources  in  the  absence  of  the  other  source, 

and  the  term  corresponds  to  the  interference  of  the  two  sources.) 

At  the  same  frequency,  the  phase  shift  ,^,  between  waves  from  two  adjacent  sources 

is  independent  of  time: 


•  When  Ii  =  I2  =  Iq,  which  is  the  case  when  the  two  sources  have  the  same 
intensity,  we  have  I  =  AIq  cos^  For  the  case  with  ^  =  0,  the  phase  shift  offsets 
the  phase  difference  due  to  the  separation  of  the  sources  and  we  have  I  =  4/o. 

•  Figures  3-2  and  3-3  show  the  interference  between  two  sources  at  four  typical 
separations  of  A,  A/2,A/4  and  A/8  with  (^0  =  0).  From  these  plots,  it  can 
be  seen  that  within  one  wavelength,  as  the  distance  of  the  sources  decrease. 
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the  intensity  becomes  less  directionally  dependent.  AU  of  the  plots  have  the 
maximum  intensity  at  0  =  0  as  they  are  always  ‘in  phase’  in  that  direction. 
For  the  case  in  which  the  separation  equals  one  wavelength,  the  sound  is  also  in 
phase  in  the  vertical  direction  and  so  there  is  a  maximum  intensity  downwards 
and  upwards.  For  a  half  wavelength  separation,  the  sources  are  exactly  out  of 
phase  and  there  is  no  response  at  0  =  90°.  For  the  case  in  which  the  separation 
is  just  an  eighth  of  the  wavelength,  the  intensity  pattern  is  very  close  to  the 
point  source  case. 


3.3  Multiple  Source  Interference 

3.3.1  Analysis  Using  Complex  Numbers 

Complex  numbers  can  simpKfy  the  analysis  of  the  interference  of  multiple  sources  (see 
Clay  and  Medwin,  1976).  Many  operations  involve  the  sums  and  differences  of  angles 
and  the  products  of  trigonometric  functions.  These  operations  are  simplified  by  using 
the  relations  between  trigonometric  functions  and  complex  exponential  functions. 


=  cos  #  +  i  sin  $  (3.11) 

cos$  =  - - - ,sin$  = - — -  (3.12) 

For  N  sources  evenly  spaced  over  a  distance  W,  the  separation  of  neighboring 
sources  is: 


d  = 


W 


N-1 


(3.13) 


and  the  pressure  fluctuation,  Ap„,  of  the  signal  from  the  nth  source,  relative  to  the 
source  at  the  distance  R,  is 


A  r-/  ,  1  nkW  sin  6., 

Apn  =  aexp[i{u;t  -  kR+ (3.14) 

where  a  is  a  constant. 
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So,  the  pressure  fluctuation  of  n  sources  is 


^■-1 


A  rv  ^  7  DM  ,inkWsm4>. 

Ap  =  aexp[%{u;t  -  kR)\  2^  exp( — - ^ — ) 


n=0 


N-l 


(3.15) 


Since  wt  —  kR  is  common  to  all  the  signals,  we  factor  it  at  the  beginning  and  then 
suppress  it  by  calculating  C  as  follows: 


^  AnkW  SID.  (j) . 

C  =  al^  exp{  _  — ) 

n=0  ^ 


N-l 

=  a  ^2  exp{inkd  SID.  4>) 

71=0 

We  can  show  that  (see  p.  46,  Clay  and  Medwin,  1976): 


smlNk(W/2)^\ 

_ r '  7  / TTr /n\  sincA  1  J 


(3.16) 


exv\ih{WI2Y^]  ‘  Nsm[k(WI2)fYi] 

The  expression  in  braces  has  an  absolute  value  of  1  and  specifies  a  phase  shift  that 
depends  on  the  choice  of  origin.  The  remaining  factor  is  known  as  the  “directional 
response”,  D,  which  is: 


D  = 


sin(^flsin7^) 

■^sin(^^sin<?i)) 


(3.17) 


When  N  is  large,  and  using 


sin[A:(W/2)(sin<^/(iV  -  1))]  A:(W/2)(sin  </>)/( W  -  1) 


(3.18) 


D  becomes 


D  = 


sin 


kW  sin0 
2 


kW  sin  4> 


(3.19) 


2 

The  latter  expression  has  the  form  (sinx)/x  which  has  a  maximum  of  one  as  x 
tends  to  zero.  This  is  identical  to  the  directional  response  of  a  continuously  distributed 
line  source  (Clay  and  Medwin,  1976). 
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If  we  substitute  S  =  kdsmcj)  [d  is  the  distance  between  two  neighboring  sources) 
into  C,  we  have 


C  —  N aex'p[ik(W 1 2)  sm  (j)] 


sin[Nk{d/2)  sin  0] 
N  sin(/^(<i/2)  sin  4>) 


rvAr  sin^ 

=  aexp[i{N  —  1)-]—  ^ 


2  sin(5/2) 

So,  we  have  the  expression  for  the  energy  intensity  of  n  equally  spaced  sources 
(denoting  Iq  =  a^) 


^  ^  ,  sin  ^ ,  ,  , 

I  =  Io*{-:^y  3.20 

sin  5 

Figures  3-4  and  3-5  show  us  the  situation  without  considering  a  free  surface.  The 
plots  have  the  following  characteristics: 

•  When  the  sources  are  in  phase  (as  dsin^  =  ihfcA),  the  intensity  reaches  its 
maximum  (/  =  N^Iq).  As  for  the  two  source  case,  when  the  sources  are  in- 
phase,  there  is  always  a  maximum  response  horizontal  to  the  vertical  array 
which  is  independent  of  the  distance  between  the  sources.  More  importantly, 
the  maximum  direction  stays  the  same  when  we  add  more  sources  to  the  array 
because  the  new  sources  are  stiU  in  phase  with  the  old  sources. 

•  Also,  we  notice  that  in  the  direction  of  maximum  intensity,  all  of  the  sources 
must  be  in  phase  with  each  other. 

•  Generally,  nodes  in  the  beam  pattern  increase  as  the  number  of  sources  in¬ 
creases.  In  theory,  at  the  position  that  6  =  k'  =  1,2, ... ,  [N  —  1),  {N  -j- 

1), . . .  ,  (under  the  condition  that  k'  ^  0)  the  intensity  reaches  its  minimum 

7  =  0. 

•  Between  the  minimum  position,  there  are  local  maximum  responses.  We  can 
get  the  local  maximum  position  and  its  intensity  theoretically  by  letting  ^  =  0. 


33 


Beam  angle  =  0  (phase  delay  =  0) 


Beam  angle  =  30  (phase  delay  =  109.8) 


Beam  angle  =  15  (phase  delay  =  56.8) 


Beam  angle  =  45  (phase  delay  =  155.3) 


Figure  3-4:  Interference  Patterns  for  a  Ten  Element  Array  without  a  Free  Surface  (I), 
plotted  as  rose  diagrams 
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Beam  angle  =  60  (phase  delay  =  1 90.2)  Beam  angle  =  75  (phase  delay  =  212.1) 


Beam  angle  =  90  (phase  delay  =  219.6)  Beam  angle  =  105  (phase  delay  =  212.1) 


Figure  3-5:  Interference  Patterns  for  Ten  Sources  without  a  Free  Surface  (II),  plotted 
as  rose  diagrams 


0  I/I  0  I/I 


Beam  angle  =  0  (phase  delay  =  0)  Beam  angle  =  15  (phase  delay  =  56.8) 


Beam  angle  =  30  (phase  delay  =  109.8)  Beam  angle  =  45  (phase  delay  =  155.3) 


Figure  3-6:  Ten  Element  Array  without  a  Free  Surface  (I),  plotted  as  intensity  versus 
angle 


Beam  angle  =  90  (phase  delay  =  219.6) 


Beam  angle  =  105  (phase  delay  =  212.1) 


Figure  3-7:  Ten  Element  Array  without  a  Free  Surface  (II),  plotted  as  intensity  versus 


This  leads  to  positions  (S  =  kd  sin  <l>)  at 


.  AT*  ^ 

tan  — —  =  N  tan  — 

2  2 


(3.21) 


ajid  the  local  maxima  of 


NH 

1  +  (iV2  -  l)sin2  f 


(3.22) 


3.3.2  Multiple  Sources  with  a  Free  Surface 

In  the  case  of  multiple  sources  at  depth  h  below  a  free  surface,  the  fluctuation  has 
two  parts,  one  is  directly  from  the  n  sources. 


../  .  T,  ...  AnkW  sin  (f>. 

Api  =  aexp[t[u}t  —  kR  +  whsin  <p)\  2_^  exp[ — — — - — ) 

n=0  iV  —  1 

and  the  other  part  corresponds  to  their  images 

.  r»/  .T.  .../\i  V" '  ^  “tTikiV  sxn  ^ 

Ap2  =  aexp[i{(jjt  —  kR  —  kh  sm  <p)\  2_,  ea;p( — — j j 


We  have  already  shown  that 


,inkW sin  <l>.  ,  ,  sin  ^ 


a  E  e;rp( 


sin(^/2) 


Similarly  we  can  show  that 


^  ^-inkW sin  <l>^  r  sin  ^ 

“  ?„  N-I  )  = 


So,  the  total  fluctuation  is 


(3.23) 


(3.24) 


(3.25) 


(3.26) 


Ap  =  Api  +  Ap2 

sin^ 
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3.3.3  Multiple  Sources  under  a  Free  Surface  with  Phase 
Delay 

By  modifying  the  phase  difference  between  two  neighboring  sources,  we  can  control 
the  outgoing  beam  direction.  Figures  3-4  and  3-5  show  the  situation  without  con¬ 
sideration  of  a  free  surface.  Figures  3-8  and  3-9  show  us  the  interference  pattern  of 
10  sources  Avith  different  outgoing  beam  angle  with  consideration  of  a  free  surface. 

Assuming  that  the  phase  difference  between  two  neighboring  sources  is  the  same, 
with  phase  advance  we  have 

A^,  =  aexp[i(wt  -kR+kh  sin  4,)]  +  i„)  (3.29) 

n=0  ^  1 
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0  I/I  0  I/I 


Beam  angle  =  30  (phase  delay  =  109.8)  Beam  angle  =  45  (phase  delay  =  155.3) 


Figure  3-8:  Interference  Patterns  for  Ten  Element  Array  with  a  Free  Surface  (I), 
plotted  as  rose  diagrams 


Beam  angle  =  90  (phase  delay  =  219.6) 


Beam  angle  =  105  (phase  delay  =  212.1) 


Figure  3-9:  Interference  Patterns  for  Ten  Element  Array  with 
plotted  as  rose  diagrams 


Free  Surface  (II), 


0  I/I  0  I/I 


Beam  angle  =  0  (phase  delay  =  0)  Beam  angle  =  IS  (phase  delay  =  S6.8) 


Beam  angle  =  30  (phase  delay  =  109.8)  Beam  angle  =  45  (phase  delay  =  155.3) 


Figure  3-10:  Ten  Element  Array  with  a  Free  Surface  (I),  plotted  as  intensity  versus 
angle 


0  m 


Beam  angle  =  90  (phase  delay  =  219.6)  Beam  angle  =  105  (phase  delay  =  212. 1) 


Figure  3-11:  Ten  Element  Array  with  a  Free  Surface  (II),  plotted  as  intensity  versus 


and  another  part  corresponding  to  their  images 


A  r*/  iT^  V  '  /  zTzk^V  sin  ^  \ 

Ap2  =  aexp[i{u)t  —  kn  —  kh  sin  (f>)\  ^  €xp( — — +  ooj 

n=0  jv  —  1 


(3.30) 


Let  b  =  +  ^0  and  b'  =  the  total  fluctuation  is  expressed  as 


Sin  2 


sin 


(3.31) 


and  the  energy  density  is 


sin  ^  sin  ^ 

/  =  j4i  +  ^2  —  2  *  \ *  cos(2kh  sin  <j>  +  kW  sin  <j>)  (3.32) 

sin|sin|- 

where 

sin=^^ 

^  *  iC-^r  (3-33) 

sin  2 

sin  ^ 

A2  =  aU  {(!^)2  (3.34) 

sm  j 

When  So  =  0,  we  get  the  case  without  phase  difierence  between  sources.  From 
equation  3.20,  we  notice  that  Ai  is  the  same  as  the  ten  element  interference  result 
without  a  free  surface.  A^  is  also  a  ten  element  interference  result,  but  it  has  a 
reversed  phase  delay  corresponding  to  the  image  of  the  vertical  array.  Also,  the  part 
—2  ♦  "  *  cos(2kh  sin  (f>  +  kW  sin  4>)  in  the  equation  3.32  can  be  considered 

sin  j  sin  y 

as  the  interference  between  the  ten  element  array  and  its  image  without  considering 
the  free  surface. 


Chapter  4 


A  Study  of  a  Ten  Element 
Vertical  Array  in  a  Depth 
Dependent  Ocean 


4.1  Program  Development  for  a  Vertical  Array 

The  weakness  of  the  reflectivity  model  is  that  it  requires  large  computation  time. 
The  original  program  has  noise  in  the  far  field,  say  beyond  30  km  offset  at  50  Hz. 
Attenuation  is  included  in  order  to  compress  the  noise,  which  alters  the  result.  The 
new  program  finds  a  way  to  improve  the  accuracy  and  also  avoids  the  use  of  attenu¬ 
ation.  In  order  to  improve  the  accuracy  of  modeling,  we  have  to  make  S'p  very  small. 
While  making  it  smaller,  we  increase  the  relative  error  in  the  truncation  process  in 
the  computation,  which  results  in  noise  in  the  final  modeling(see  section  2.3).  We  get 
rid  of  the  noise  by  successfully  compressing  the  truncation  error  in  the  far  field.  The 
new  program  is  optimized  both  in  speed  and  structure  of  the  program.  Also,  a  new 
plotting  program  was  written  in  matlab  4.0  to  plot  the  amplitude  versus  distance. 


4.2  New  Program  Flow  Chart 

The  ocean  medium  can  be  assumed  to  be  stratified.  The  sound  speed  in  the  water 
column  is  approximated  by  thin  homogeneous  layers  with  small  discontinuous  jumps 
in  velocity  across  interfaces.  The  new  program  has  the  following  flow  chart: 

Sources  Loop 

get  the  neH  source  depth; 

determine  the  layer  number  that  the  source  belong  to; 
change  the  model  so  that  the  source  is  on  the  top  the  layer; 

get  the  receivers  depths; 

determine  the  layer  number  that  the  receivers  belong  to; 
change  the  layer  so  that  the  receivers  are  on  the  top  the  layer; 

Frequency  Loop 
Layer  Loop 
p  loop 

Iteration  equations  to  confute  the  doimweu:d(eilso  upweurd) 
looking  reflection  emd  tremsmission  coefficient  matrices . 
end-Cp  loop} 

If  source  layer  then  {vrite  down  its  reflection 
and  transmission  coefficient  matrices.  } 

If  receiver  layer  then  {write  down  its  reflection 
and  tremsmission  coefficient  matrices.  }■ 
end{Layer  Loop} 

Receiver  Loop 
X  loop 
p  loop 

integrate  over  p  to  transform  p  to  x 
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end-Cp  loop} 
end  {x  loop} 
end-CReceiver  Loop} 

multiply  the  phase  delay  factor  to  the  source  response 
end-CFrequency  Loop} 

sum  the  new  soxirce  response  to  total  update  sources  response 
end-CSources  Loop} 

It  is  easy  to  get  the  vertical  array  response.  We  simply  edit  a  file  which  has  sources 
and  receiver  information.  Generally,  it  works  as  follows: 

•  Get  the  new  source  depth,  which  is  always  deeper  than  the  previous  one.  The 
distance  between  sources  stays  the  same  and  so  the  new  source  is  always  the 

•  same  distance  deeper  than  the  previous  one. 

•  Determine  the  layer  numbers  for  the  new  sources  and  receivers.  The  program 
requires  that  the  sources  and  receivers  should  be  on  the  top  of  a  layer.  It  changes 
the  thickness  of  neighboring  layers  so  that  adding  sources  will  not  affect  the 
other  parts  of  the  layer  model. 

•  Get  the  source  response. 

•  Multiply  the  response  by  the  phase  delay  factor  to  get  the  response  relative  to 
the  top  source. 

•  Sum  up  all  of  the  responses  to  get  the  multiple  source  response. 

4.3  Reflectivity  Modeling 

4.3.1  Speed  of  Sound  and  Ray  Paths  in  the  Ocean 

The  sound  travel  path  is  associated  with  the  dependence  of  the  sound  speed  on  depth. 
Simplifying  the  results  of  the  Naval  Research  Laboratory  studies  (Medwin,  1975),  the 


47 


Figure  4-1:  Ray  diagram  for  typical  Atlantic  Ocean  sound  channel,  depicting  chan¬ 
neled  rays  and  refracted-surface-reflected(RSR)  rays;  sound  speed  profile  is  at  the 
right.  The  angles  are  grazing  angles  at  the  axis  of  the  sound  channel.  (Ewing  and 
Worzel,  1948.) 


dependence  of  sound  speed  on  temperature,  salinity,  and  pressure  within  1  km  of  the 
surface  is  found  to  be  approximately 


c  =  1449.2-|-4.6T-0.055r2-t-0.00029r®-f(1.34-0.010r)(S’-35)-|-1.58xl0-®Pa  (4.1) 


where  c  =  sound  speed  (m/s),  T  =  temperature  (°C7),  S  =  salinity  (ppt),  Pa  = 
gauge  pressure  due  to  water  column,  (N/m^). 

Ignoring  compressibility  Pa  =  paQz  and  using  pA  (1  +  5  X  10  ®)  kg/m^,  g  — 
9.8m/s*,  and  z  =  depth(m),  the  pressure  term  can  be  obtained. 

When  the  temperature  has  a  large  decrease  with  increa.sing  depth,  the  temperature 
effect  overrides  the  pressure  effect  and  the  sound  speed  gradient  is  negative,  otherwise, 
the  sound  gradient  is  positive. 

Tliis  vertical  stratification  of  the  deep  ocean  and  the  resulting  ray  paths  depends 
on  location,  season,  and  time  of  day.  The  Figure  4-1  shows  the  ray  traces  computed 
by  Ewing  and  Worzel  (1948)  as  well  as  the  sound  profile.  The  steepest  ray  shown  just 
grazes  the  bottom.  The  angles  are  measured  relative  to  the  horizontal  in  long-range 
ray  tracing  and  are  called  “grazing  angles”. 
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SOUND  SPEED  (m/sec) 


Figure  4-2:  The  profile  of  sound  speed  used  in  the  modeling 


4.3.2  Model  Description 

In  Figure  4-2,  we  show  the  ten  element  vertical  array  position  and  the  ocean  velocity 
profile.  A  total  of  172  layers  are  used  to  approximate  the  smooth  gradient  zones.  We 
put  the  sources  in  such  a  way  that  all  of  the  sources  are  in  different  layers.  In  this 
chapter  we  present  results  at  50  Hz  (compared  to  250  Hz  in  the  last  chapter).  At  this 
frequency  the  refiectivity  solutions  are  faster  and  have  less  noise. 

4.3.3  Homogeneous  Result 

The  results  for  a  homogeneous  ocean  are  shown  in  figures  from  Figure  4-3  to  Figure  4- 
9.  They  show  that: 

•  The  ‘deeper  the  receivers  are,  the  closer  the  modeling  result  agrees  with  the 
■  theoretical  result.  We  compare  the  modeling  residts  at  250  meters  below  the 

top  source  of  the  vertical  array  with  the  theoretical  results,  and  find  that  they 
are  similax  in  general  characteristics  but  there  are  differences  in  the  local  chax- 
acteristics.  For  receivers  at  4000  m  depth  the  modeling  results  agree  well  with 
theoretical  predictions.  The  reason  for  the  difference  in  the  receiver  responses  at 
250  m  and  4000  m  can  be  deduced  from  the  conditions  we  put  on  the  theoretical 
work.  We  assumed  that  the  receivers  were  an  infinite  distance  from  the  vertical 
array.  The  receivers  at  4000  m  depth  are  much  closer  to  that  assumption  than 
those  at  250  m.  We  conclude  that  the  theoretical  results  hold  pretty  weU  in  the 
far  field  but  axe  less  accurate  in  the  near  field. 

•  Comparing  the  reflectivity  results  with  the  theoretical  results,  we  must  consider 
the  geometric  effect  on  the  responses.  The  theoretical  results  assume  that  all  of 
the  receivers  axe  at  the  same  distance  from  the  source,  while  in  our  models  all  of 
the  receivers  are  laid  out  horizontally.  So,  because  of  the  geometric  spreading 
effect,  the  responses  of  receivers  on  the  sea  floor  will  be  the  product  of  the 
theoretical  results  and  a  geometric  spreading  factor  which  is  cos(angleo)-  The 
angleo  is  the  angle  between  vertical  axis  and  the  receiver  angle.  Because  the 
distance  from  the  source  to  the  receiver  is  d  =  ho/cos{angleo),  the  response 
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decreases  with  the  distance  from  the  receiver  to  the  source.  The  response  at 
distance  has  to  be  multiplied  by  the  factor  of  cos(angleo)  to  get  the  real  response 
on  the  horizontal  receivers. 

•  For  results  beyond  30°  dip  angle,  the  energy  focuses  mainly  between  zero  offset 
and  10  km  offset.  The  far  offset  response  is  barely  noticeable.  As  the  beam  be¬ 
comes  directed  downwards,  there  are  less  variations  in  the  near  offset  response. 
The  sidelobe  has  less  effect. 

•  We  do  not  notice  much  difference  in  the  results  between  the  60°,  75°  and  90° 
beams  at  4000Tn,  except  that  amplitude  decreases  as  the  beam  angle  increases. 
The  reason  that  the  maximum  mainly  remains  at  1-3  km  distance  is  that  in  the 
far  field,  the  response  results  from  sidelobes  of  the  interference  pattern. 

•  The  vertical  array  response  can  be  divided  into  three  components:  the  vertical 
array  response  without  free  surface,  the  image  of  the  vertical  array  response  in 
the  homogeneous  media,  and  the  interference  of  the  vertical  array  and  its  image 
array.  The  complex  interference  between  the  array  and  its  image  introduces 
the  ‘finger’  feature  into  the  response  on  the  seafloor.  The  results  show  that  the 
image  of  the  array  response  gives  a  better  fit  to  the  vertical  array  response  than 
the  response  from  the  array  itself.  Also,  the  theoretical  vertical  array  response 
shows  stronger  ‘fingering’  than  the  modeling  results. 

4.3.4  Gradient  Ocean  Result 

The  results  for  the  ocean  with  gradient  velocity  profile  are  shown  as  follows: 

•  The  gradient  effect  can  be  seen  from  Figure  4-11  to  Figure  4-16.  Unlike  the  ho¬ 
mogeneous  situation,  the  minimum  response  never  goes  to  zero.  This  is  because 
the  sources  are  at  depths  with  different  velocities  and  all  of  the  waves  reaching 
the  receivers  have  complex  paths.  In  this  case  the  assumptions  for  the  theory 
of  a  homogeneous  ocean  can  not  be  applied. 
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•  The  gradient  has  the  effect  of  increasing  the  amplitudes  between  25  and  30  km 
by  about  6  dB  to  the  0®,  15°,  30°  beams.  For  the  rest  of  the  offsets,  the  modeling 
results  for  receivers  at  a  depth  of  40  007n  show  a  good  fit  to  the  homogeneous 
theoretical  results.  This  indicates  that  for  our  arrangement  of  a  vertical  array, 
in  a  general  ocean  velocity  gradient,  the  beam  is  affected  at  some  ranges  and 
in  our  case  there  is  a  great  increase  in  amplitude  between  20  and  30  km  range, 
especially  for  0°,  15°  and  30°  beams.  For  the  rest  of  the  offsets,  the  beam  does 
not  change  its  characteristics  dramatically.  We  conclude  that  the  gradient  in 
velocity  profile  has  a  big  impact  on  the  beam  pattern. 

•  We  found  our  results  are  quite  different  from  the  work  done  by  Ozliier  using  the 
ray  tracing  technique,  even  considering  that  we  used  different  oceanic  velocity 
models.  First,  there  is  no  clear  indication  of  the  interference  between  sources 
and  its  image  on  her  results  which  should  change  the  beam  pattern  dramatically. 
Second,  her  resiilts  found  no  indication  that  the  oceanic  velocity  profile  has 
effects  on  the  beam.  Third,  our  modeling  results  show  interference  even  beyond 
50  km  offset.  Her  restdts  cannot  be  trusted  after  28  km  offset.  Fourth,  our 
results  show  patterns  due  to  the  complicated  interference  between  sources  and 
their  images,  her  results  show  a  pretty  simple  interference  pattern.  Also,  our 
model  uses  50  Hz  sources  and  her  model  used  250  Hz  sources  which  shoTild  have 
much  strong  interference  patterns.  It  might  be  helpful  if  we  can  redo  her  work 
applying  our  program  to  the  same  model.  This  can  be  part  of  future  work. 

•  It  would  be  interesting  to  investigate  the  vertical  array  resiilts  of  a  pulse.  Also, 
studjdng  the  wave  path  using  ray  tracing  techniques  might  give  us  more  insight 
into  the  vertical  array  response  for  the  gradient  ocean.  It  is  believed  that  several 
velocity  models  might  need  further  study. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-3:  The  homogeneous  ocean  response  for  a  10  element  vertical  array  under  a 
free  surface  with  a  0  degree  beam  angle.  The  agreement  is  good  out  to  40  km.  The 
solid  line  is  for  the  reflectivity  modeling,  the  dash-dot  line  is  for  theoretical  results, 
the  dash  line  is  the  theoretical  result  for  the  ten  element  array  and  also  its  image 
(they  are  the  the  same).  Note  that  all  curves  have  been  normalized  to  a  maximum 
value  of  one. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-4:  The  homogeneous  ocean  response  for  a  10  element  vertical  array  under  a 
free  surface  with  a  15  degree  beam  angle.  The  agreement  is  good  out  to  40  km.  The 
solid  line  is  for  the  reflectivity  modeling,  the  dash-dot  line  is  for  theoretical  results, 
the  dash  line  and  dots  line  are  the  theoretical  residts  for  the  ten  element  array  and  its 
image  respectively.  Note  that  all  curves  have  been  normalized  to  a  maximum  value 
of  one. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-5:  The  homogeneous  ocean  response  for  a  10  element  vertical  array  under  a 
free  surface  with  a  30  degree  beam  angle.  The  agreement  is  good  out  to  40  km.  The 
solid  line  is  for  the  reflectivity  modeling,  the  dash-dot  line  is  for  theoretical  results, 
the  dash  line  and  dots  line  are  the  theoretical  results  for  the  ten  element  array  and  its 
image  respectively.  Note  that  all  curves  have  been  normalized  to  a  maximum  value 
of  one. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-6:  The  homogeneous  ocean  response  for  a  10  element  vertical  array  under  a 
free  surface  with  a  45  degree  beam  angle.  The  agreement  is  good  out  to  40  km.  The 
solid  line  is  for  the  reflectivity  modeling,  the  dash-dot  line  is  for  theoretical  results, 
the  dash  line  and  dots  line  are  the  theoretical  results  for  the  ten  element  array  and  its 
image  respectively.  Note  that  all  curves  have  been  normalized  to  a  maximum  value 
of  one. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-7:  The  homogeneous  ocean  response  for  a  10  element  vertical  array  under  a 
free  surface  with  a  60  degree  beam  angle.  The  agreement  is  good  out  to  40  km.  The 
solid  line  is  for  the  reflectivity  modeling,  the  dash-dot  line  is  for  theoretical  results, 
the  dash  line  and  dots  line  are  the  theoretical  results  for  the  ten  element  array  and  its 
image  respectively.  Note  that  all  curves  have  been  normalized  to  a  maximum  value 
of  one. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-8:  The  homogeneous  ocean  response  for  a  10  element  vertical  array  under  a 
free  surface  with  a  75  degree  beam  angle.  The  agreement  is  good  out  to  40  km.  The 
solid  line  is  for  the  reflectivity  modeling,  the  dash-dot  line  is  for  theoretical  results, 
the  dash  line  and  dots  line  are  the  theoretical  results  for  the  ten  element  array  and  its 
image  respectively.  Note  that  all  curves  have  been  normalized  to  a  maximum  value 
of  one. 


58 


Response  at  depth  4000  m  below  the  top  source 


Figure  4-9:  The  homogeneous  ocean  response  for  a  10  element  vertical  array  under  a 
free  surface  with  a  90  degree  beam  angle.  The  agreement  is  good  out  to  40  km.  The 
solid  line  is  for  the  reflectivity  modeling,  the  dash-dot  line  is  for  theoretical  results, 
the  dash  line  and  dots  line  are  the  theoretical  results  for  the  ten  element  array  and  its 
image  respectively.  Note  that  aU  curves  have  been  normalized  to  a  maximum  value 
of  one. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-10:  The  gradient  ocean  response  for  a  10  element  vertical  array  (50  Hz) 
with  a  0  degree  beam  angle  (solid  line).  The  gradient  has  the  effect  of  increasing  the 
amplitudes  between  25  and  30  km  by  about  6  dB.  Because  of  geometrical  spreading 
the  largest  response  at  the  seafloor  occurs  at  8  km  from  side  lobes.  The  dashed  line  is 
the  theoretical  result  for  the  ten  element  array  and  is  given  as  a  reference.  Note  that 
all  curves  have  been  normalized  to  a  majdmum  value  of  one  and  it  would  be  useful 
to  compare  reflectivity  for  both  homogeneous  and  gradient  models. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-11:  The  gradient  ocean  response  for  a  10  element  vertical  array  (50  Hz) 
with  a  15  degree  beam  angle  (solid  line).  The  gradient  has  the  effect  of  increasing 
the  amplitudes  between  25  and  30  km  by  about  6  dB.  For  a  homogeneous  ocean  a  15 
degree  beam  would  intersect  the  seafloor  at  15  km.  However  because  of  geometrical 
spreading  the  largest  response  at  the  seafloor  occurs  at  8  km  from  side  lobes.  The 
dashed  line  is  the  same  as  dash-dot  line  in  the  homogeneous  ocean  response  plot  and 
is  given  as  a  reference.  Note  that  all  curves  have  been  normalized  to  a  maximum 
value  of  one  and  it  would  be  useful  to  compare  reflectivity  for  both  homogeneous  and 
gradient  models. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-12:  The  gradient  ocean  response  for  a  10  element  vertical  array  (50  Hz) 
with  a  30  degree  beam  angle  (solid  line).  The  gradient  has  the  effect  of  increasing 
the  amplitudes  between  25  and  30  km  by  about  6  dB.  The  dashed  line  is  the  same 
as  dash-dot  line  in  the  homogeneous  ocean  response  plot  and  is  given  as  a  reference. 
Note  that  all  curves  have  been  normalized  to  a  maodmum  value  of  one  and  it  would 
be  useful  to  compare  reflectivity  for  both  homogeneous  and  gradient  models. 
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Figure  4-13:  The  gradient  ocean  response  for  a  10  element  vertical  array  (50  Hz) 
with  a  45  degree  beam  angle  (solid  line).  The  gradient  has  the  effect  of  increasing 
the  amplitudes  between  25  and  30  km  by  about  6  dB.  The  dashed  line  is  the  same 
as  dash-dot  line  in  the  homogeneous  ocean  response  plot  and  is  given  as  a  reference. 
Note  that  all  curves  have  been  normalized  to  a  maximum  value  of  one  and  it  would 
be  useful  to  compare  reflectivity  for  both  homogeneous  and  gradient  models. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-14:  The  gradient  ocean  response  for  a  10  element  vertical  array  (50  Hz) 
with  a  60  degree  beam  angle  (solid  line).  The  gradient  has  the  effect  of  increasing 
the  amplitudes  between  25  and  30  km  by  about  6  dB.  The  dashed  line  is  the  same 
as  dash-dot  line  in  the  homogeneous  ocean  response  plot  and  is  given  as  a  reference. 
Note  that  all  curves  have  been  normalized  to  a  maximum  value  of  one  and  it  would 
be  useful  to  compare  reflectivity  for  both  homogeneous  and  gradient  models. 


64 


Response  at  depth  4000  m  below  the  top  source 


Figure  4-15:  The  gradient  ocean  response  for  a  10  element  vertical  array  (50  Hz) 
with  a  75  degree  beam  angle  (solid  line).  The  gradient  has  the  effect  of  increasing 
the  amplitudes  between  25  and  30  km  by  about  6  dB.  The  dashed  line  is  the  same 
as  dash-dot  line  in  the  homogeneous  ocean  response  plot  and  is  given  as  a  reference. 
Note  that  aU  curves  have  been  normalized  to  a  maximum  value  of  one  and  it  would 
be  useful  to  compare  reflectivity  for  both  homogeneous  and  gradient  models. 
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Response  at  depth  4000  m  below  the  top  source 


Figure  4-16:  The  gradient  ocean  response  for  a  10  element  vertical  axray  (50  Hz) 
with  a  90  degree  beam  angle  (solid  line).  The  gradient  has  the  effect  of  increasing 
the  amplitudes  between  25  and  30  km  by  about  6  dB,  The  dashed  line  is  the  same 
as  dash- dot  line  in  the  homogeneous  ocean  response  plot  and  is  given  as  a  reference. 
Note  that  aU  curves  have  been  normalized  to  a  maximum  value  of  one  and  it  would 
be  useful  to  compare  reflectivity  for  both  homogeneous  and  gradient  models. 
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Chapter  5 

Comparison  of  Modified 
Refiectivity  Results  with 
Benchmark  Solutions 


5.1  Test  Problem  Solutions 


In  order  to  test  our  reflectivity  code,  we  compare  solutions  of  two  of  seven  test  prob¬ 
lems  discussed  at  the  PE  Workshop  II  held  in  Slidell,  LA  on  6-10  May  1991  (Jensen, 
1993). 


5.2  Test  Case  1 

A  point  source  placed  near  a  reflecting  boundary  in  a  homogeneous  medium  gives 
rise  to  the  well-known  Lloyd-mirror  interference  pattern  for  which  an  exact  field  so¬ 
lution  is  available.  The  Lloyd-mirror  pattern  is  ideal  for  checking  both  the  angular 
distribution  of  energy  associated  with  a  given  starting  field  and  the  high-angle  capa¬ 
bility  of  parabolic  methods.  It  is  also  a  simple  test  for  checking  the  normalization  of 
reflectivity  output  to  transmission  loss. 

In  Test  Case  1  the  fluid  halfspace  has  a  constant  speed  of  1500  m/s.  The  source 
frequency  is  40  Hz  and  the  source/receiver  depths  below  the  free  surface  are  350m  and 
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3990m.  The  solid  line  in  Fig.  5-1  shows  the  image  reference  solution  which  is  known 
to  be  exact  reference.  The  reflectivity  result  is  in  perfect  agreement  with  the  reference 
solution  out  to  10  km  range  at  40  Hz.  Numerically  stable  results  were  obtained  with 
Sp  =  0.00005,  starting  p  =  0.0,  ending  p  =  1.0.  In  order  to  get  this  agreement  a  DC 
offset  of  6  dB  was  applied  to  the  reflectivity  result. 

5.3  Test  Case  7 

During  PE  Workshop  II,  test  case  7  was  introduced  which  reveals  serious  deficiencies 
in  some  PE  approximations.  The  test  problem  deals  with  long-range  propagation  in  a 
surface  duct.  The  source  frequency  is  80  Hz,  which  results  in  a  leaky  (virtual)  mode, 
which  continuously  sheds  energy  into  the  lower  medium,  as  indicated  by  the  dashed 
arrows  in  Figure  7.1  of  Jensen,  1993. 

The  test  problem  is  critical  to  test  that  a  small  phase  error  for  the  refracted  leaky 
path  can  result  in  large  changes  in  sound  level  in  the  duct  beyond  the  first  convergence 
zone.  Thus,  in  the  case  where  the  two  arrivals  (the  ducted  and  the  leaky  arrival)  have 
similar  amplitudes,  we  may  observe  level  changes  of  several  tens  of  dB’s  between 
constructive  and  destructive  path  interference. 

The  reference  solution  to  this  problem  for  a  source  at  25  m  and  a  receiver  at  100 
m  was  generated  by  the  SNAP  normal- mode  code.  The  SNAP  reference  solution  is 
shown  as  the  solid  line  in  the  Fig.  5-2.  The  reflectivity  results  are  shown  as  a  dashed 
line.  Numerically  stable  results  were  obtained  with  a,  Sp  =  0.00005,  starting  p  =  0.0, 
ending  p  =  1.0.  The  number  of  layers  was  760  which  makes  a  layer  thickness  of  about 
A/4  .  The  results  are  poorer  at  coarser  thicknesses.  When  starting  receiver  range  = 
0.5  km,  ending  receiver  range  =  150.5  km,  receiver  offset  =  0.81  km.  The  run  time 
on  an  IBM  RS-6000  is  22  minutes  and  10  seconds  for  one  frequency  (80  Hz). 

The  results  show  that  the  our  reflectivity  code  can  handle  this  problem  accurately 
for  80  Hz  to  range  beyond  100  km.  No  offset  was  applied  to  the  reflectivity  output 
to  get  this  agreement. 
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Amplitude  (dB) 


Figure  5-1:  Comparison  of  image  reference  solution  (solid  line)  with  reflect] 
suit  (dotted  line) 
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SNAP/Reflectivity  Method  Solution 


Figure  5-2:  Compaxison  of  normal-mode  reference  solution  (SNAP  with  solid  line) 
with  reflectivity  result(dotted  line) 
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5.4  Discussion 


In  order  to  get  correct  response  in  the  far  field,  I  modified  the  original  program  to  get 
the  acurate  response  in  far  field.  There  is  the  several  improvement  in  the  reflectivity 
code. 

•  I  change  the  single  precision  program  to  double  precision  which  increases  the 
the  precision  of  the  program  and  thus  reduces  the  error  of  reflectivity  result. 

•  The  new  program  is  much  simpler  which  help  to  reduce  unnecessary  numerically 
error.  The  new  program  can  interpolate  the  velocity  profile  in  a  easy  way  so 
that  it  is  easy  to  decrease  the  layer  thickness  until  the  result  convergence.  In 
our  practice,  we  found  the  reflectivity  results  converge  when  the  layer  thickness 
is  about  a  quarter  of  the  wavelength. 

•  It  is  essential  to  get  the  correct  result  by  introducing  an  attenuation  factor 
in  the  frequency,  which  the  program  uses  to  get  around  the  singularity  pole 
on  the  real  axis  by  deforming  the  integration  path.  Without  attenuation,  we 
cannot  get  an  convergence  result  for  the  PE  test  case  7.  The  attenuation  factor 
is  introduced  by  making  the  frequency  /o  complex,  which  gets  new  frequency 
/  =  (fo,  decay l{2  *  ir)),  the  decay  is  the  attenuation  factor.  In  our  calculation, 

I  found  we  can  get  quick  convergent  result  by  letting  decay  =  Zop(l.l)/4. 

•  The  new  program  correct  attenuation  effect  on  the  far  field  response  by  mul¬ 
tiplying  a  correction  factor,  where  x  is  the  horizontal  offset.  I  found 

the  factor  from  the  equation  2.33,  which  brings  a  decay  factor  after  we  apply 
complex  frequency  by  bringing  in  attenuation. 

•  The  original  program  does  not  get  amplitude  relative  to  Im  from  the  source 
which  implies  that  we  might  use  a  “magic  number”.  In  PE  Workshop  test  7, 
we  get  the  results  exactly  relative  to  Im  from  the  source.  Unfortunately,  for 
the  Test  case  1,  we  get  to  twice  big  responses  applying  our  modified  program 
In  our  future  study,  we  will  try  to  calibrate  the  source  level  to  get  amphtude 
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relative  to  Im.  I  believe  we  might  need  a  factor  of  frequency  difference  in  the 
program,  because  in  the  PE  Workshop  test  7  the  source  frequency  is  40  Hz 
instead  of  80  Hz  used  in  PE  Workshop  test  7. 

* 

•  Appl3dng  the  improved  program,  I  get  the  very  far  offset  response  (up  to  150  km) 
of  a  source  in  the  deep  ocean  at  frequency  of  250  Hz,  while  original  program 
cannot  get  the  response  without  significant  noise  beyond  50  km  for  50  Hz. 
Figure  5-3  shows  that  response  of  the  250  Hz  source  on  the  250  m  and  4000  m 
receivers  below  the  source  using  the  PE  Workshop  test  7  model  except  that  it 
uses  250  Hz  source  instead  of  80  Hz.  It  is  surprising  that  it  do  not  have  any 
noise  as  far  as  150  km  offset.  In  the  result,  we  notice  that  the  convergence 
zones  are  at  about  30  km,  75  km,  and  140  km  offset  which  confirms  Figure  7.1 
of  Jensen  1993.  It  needs  to  mention  that  the  convergence  result  is  gotten  when 
I  apply  the  layer  thickness  as  fine  as  a  quarter  of  wavelength  which  is  finer  than 
that  of  80  Hz  in  the  PE  Workshop  test  7. 
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Range  (km) 


Figure  5-3:  Reflectivity  Result  of  250  Hz  Source  for  Test  7  Velocity  Model 


Chapter  6 


Conclusions 


A  theoretical  study  has  been  carried  out  for  a  vertical  array  with  the  capability  of 
steering  the  beam  angle  in  a  homogeneous  ocean  and  in  a  homogeneous  ocean  with  a 
free  surface.  A  modified  reflectivity  technique  has  been  implemented  with  the  ability 
to  calculate  the  total  deep  ocean  response  of  a  vertical  array.  Also,  the  modified  code 
can  get  results  with  less  noise,  especially  in  the  far  field,  which  enables  us  to  get 
the  response  as  far  as  55  km  offset  for  50  Hz  sources  with  little  noise.  This  enables 
us  to  study  the  response  for  different  beam  angles  at  long  offsets.  AU  of  the  results 
axe  compared  with  theoretical  results.  For  the  verticeJ  array  with  different  beam 
angles,  the  reflectivity  results  for  a  deep  ocean  response  fit  well  with  our  theory  for  a 
homogeneous  ocean  response  which  assumes  that  the  receivers  are  at  infinite  distance 
from  the  array. 

The  vertical  array  response  can  be  divided  into  three  components:  the  vertical 
array  response  without  free  surface,  the  image  of  the  vertical  array  response  in  the 
homogeneous  media,  and  the  interference  of  the  vertical  array  and  its  image  array. 
The  complex  interference  between  the  array  and  its  image  introduces  the  ‘finger’ 
feature  into  the  response  on  the  seafloor.  Also,  geometric  spreading  has  a  big  effect 
on  the  deep  ocean  response.  As  an  example,  the  maximum  response  of  the  vertical 
array  with  zero  beam  angle  is  vertically  below  the  array,  where  it  is  closest  to  the 
vertical  array  even  though  the  response  is  from  the  sidelobe  of  the  beam. 

The  gradient  has  a  significant  effect  on  deep  ocean  responses  of  the  vertical  array. 
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From  our  study  of  the  vertical  array  with  0°,  15°  and  30°  beam  angle,  we  find  the 
response  amplitudes  between  25  and  30  km  increased  by  about  6  dB,  the  energy  is 
steered  to  that  range  due  to  the  velocity  gradient.  This  is  different  from  the  conclusion 
of  Ozliier  (  1992  )  who  applied  ray  tracing  techniques. 

It  would  be  interesting  to  investigate  the  vertical  array  results  of  a  pulse.  Also, 
studying  ray  paths  of  each  source  in  the  array  using  ray  tracing  techniques  might 
give  us  more  insight  into  the  vertical  array  response  for  the  gradient  ocean.  Also,  it 
is  believed  that  different  velocity  gradient  profiles  and  frequencies  need  future  study. 
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